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Abstract 

We discuss the recent application to strongly disordered systems of the Critical 
Minimum Energy Subspace (CrMES) method, used to limit the energy subspace 
of the Wang-Landau sampling. We compare with our results on the 3D Random 
Field Ising Model obtained by a multi-range Wang-Landau simulation in the whole 
energy range. We point out at some problems that may arise when applying the 
CrMES scheme to models having a complex free energy landscape. 



1 Introduction 



The study of phase transitions via Monte Carlo simulations has recently re- 
gained new interest due to the introduction of the so called "extended ensemble 
methods", like multicanonical simulations pj, broad histogram thecniques 
entropic sampling [3], Wang-Landau method [1], among others. These new 
thechniques increased significantly the accuracy of the Monte Carlo studies 
in the cases where large free energy barriers separate different wells in the 
free energy landscape. This step forward in the developement of the Monte 
Carlo method is analog to the already largely used clustering methods [5J [6J, 
which help to overcome the critical slowing down observed at second order 
transitions. 
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Following the arrival of these new methods an avalanche of thechniques were 
developped to further improve accuracy and performance of the different al- 
gorithms. Each one of these has been validated through the application to 
well known systems, typically those that can be solved exactly like the 2d 
Ising model or those whose numerical or aproximated analytical results are 
out of discussion like the Potts model for q = 5 (to test the case of first order 
transitions). 

Recently, Malakis, Peratzakis and Fytas [7j introduced an interesting approx- 
imation method to allow for an easy extension of a Wang-Landau (WL) [1] 
study to large systems. It is based on the fact that in the thermal averages 
sums at a given temperature T, only some terms give a relevant contribu- 
tion: those whose energies correspond to the interval around the maximum of 
the probability distribution of the energy at the corresponding temperature, 
Pt{E). Hence an algorithmic procedure to determine a restricted interval AE 
centered in the energy E of the maximum of Pt{E), called minimum energy 
subspace (MES) is proposed. In this way the WL sampling needs to be per- 
formed only in this restricted energy interval, improving the efficiency of the 
algorithm and diminishing the errors introduced when joining together the 
different parts of the density of states in a multirange simulation [S]. 

This algorithm is based on the equivalence of thermodynamical ensembles and 
on the central limit theorem: the energy probability distribution at a given 
temperature T approaches a Gaussian. Moreover, for a continuous transition, 
it is supposed that this remains the case even at Tc{L), the critical tempera- 
ture of a sample whose linear size is L. Then one expects that the width of the 
critical MES (CrMES) will be of the same order that the standard deviation 
of the energy: 



where volume of the system is given by = L'^, with d the space dimension, 
T is the temperature and C the specific heat. 

The CrMES is iteratively built, starting from the central value E and extend- 
ing the interval on both sides of it, until the difference between the specific 
heat calculated using the whole energy interval (or the exact one, if known) 
and the one calculated using the restricted iterated interval becomes less than 
a given error. Assuming that one imposes the same level of error for all the 
sizes, a lattice size dependence remains: the center of the CrMES E and its 
boundaries are functions of L. At the critical temperature, using the scaling 
law for the specific heat one gets: 
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In this way, performing a WL simulation over the whole energy range for 
small lattices and determining the CrMES for this small size will be enough 
to extrapolate the CrMES where the WL simulation has to be carried out for 
larger lattices. 

This method has been succesfuUy applied to the study of the pure 2D and 3D 
Ising model. 

In this work we are interested in the application of this technique to highly 
disorder systems, as has been proposed in [9] for the 3D Random Field Ising 
Model (3D-RFIM) with bimodal distribution of the random fields. The aim of 
this work is to show, by the study of the 3D-RFIM, why the generalization of 
the CrMES method to highly disordered systems is not straightforward. 

The article is organised as follows in Sec. II we describe the method used in 
this work, in section III we present our results and in Sec. IV we discuss how 
these results point out to the aspects of the proposed method that have to be 
handled with extreme care when applying it to the study of highly disordered 
systems. 



2 Description of model and the method 



The 3D-RFIM is one of the simplest disordered systems; its hamiltonian is 
given by: 

H = - J SiSj - hiSi (3) 

<ij> i 

where hi are the local fields of intensity which we assume to be distributed 
as follows: 

p[hi) = ]^[5{K-hQ) + 5{K + h^)] (4) 

The nature of the transition of this model is a subject of controversy since 
long ago. In [10] we have studied it using the multi-range version of the WL 
simulation on the whole relevant energy range. Our results for high values of 
the random field intensity show strong first order properties. The same result 
had also been found in [1T| using a completly different calculation method, 
namely the canonical Histogram Monte Carlo and similar features have also 
been found in the case of a gaussian distribution of random fields [12] 

To analyse the method proposed in [9] we start from the density of states 
(DOS) calculated using a multi-range WL simulation on the whole energy 
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space [To] and we recalculate thermal averages limiting our DOS data only to 
the energy subintervals corresponding to those indicated in [9]. 

With this procedure we are able to point out to the different aspects that 
should be handled with special care in order to apply this technique to highly 
disordered systems. 



3 Discussion of results 



The fact that the nature of the transition of the 3D RFIM is still a matter of 
controversy calls for a special attention on the validity of the hypotesis of the 
methods used to study this model. 

Let's discuss first the relationship between the nature of the transition and the 
CrMES technique. The method introduced in [7] is based on the assumption 
that one deals with a second order transition. This is required for Equation [2] 
and the gaussian-like shape of the Pt* (E) hypothesis to be valid. 

Our WL simulation on the whole energy range shows, in agreement with other 
works [1] [12], that for strong disorder, and big enough lattices, the Cl(T) 
curves may present a multiplicity of peaks. In order to estimate the transition 
temperature T*, we observe that one of the maxima of Cl{T), located at T*, 
is associated to a Pt*{E) curve which shows two well separated peaks of equal 
height. In general, this corresponds to the highest maximum of Cl{T). 

The double peaked Pt*{E) curves are an indication of the first order character 
of the transition, at the transition temperature T* . In very rare cases it has 
been found that the double peak is a finite size effect [13] . Recently it has been 
reported that this could be the case for the 3D- RFIM [13], though further 
calculations involving sizes beyond L = 32 are needed. 

In any case the energy probability distribution at the transition temperature 
T* for a given finite L, cannot be approached by a single gaussian. In fact its 
central value E corresponds then to the minimum of Pt*{E). Moreover, in 
the case of first order transitions Eqs.[T]and[2]are not valid, and the probability 
distribution may be approached by a double gaussian [15]. 

A modified version of the CrMES method adapted to a system undergoing a 
first order transition has been proposed in [16] . It is worthwhile noticing that 
this work deals with a system presenting geometrically frustrated interactions 
but no quenched disorder. 
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Fig. 1. L = 10, /lo = 2.1 (a) Specific heat vs. temperature, (b) Probability density 
for T=3, corresponding to the maximum of (a). 

It is clear that when studying systems where the free energy landscape has a 
complex structure due to the existence of a quenched disorder, the possibility 
of a crossover to a first order transition must be kept in mind. 

Moreover finite size effects and sample to sample fluctuations are sensitive 
points studying highly disordered systems. 

As it is often the case, the signatures of a flrst order transition are observed 
only for large enough sizes. For smaller sizes the transition looks continuous, 
as the correlation lenght may easily reach the size of the system. 

Figured] illustrates this situation for a fleld Hq = 2.1. In (a) one can see the spe- 
ciflc heat curve for a sample of linear size L = 10 and in (b) the corresponding 
Pt{E). There is only one peak in Cl{T) at T* and the corresponding proba- 
bility distribution shows a single "gaussian-like" peak. The appealing idea of 
the CrMES technique should be considered step by step in this case. In [7] 
it is proposed to carry WL simulations in the whole energy range on lattices 
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Fig. 2. /iQ = 2.1 Specific heat vs. temperature, comparison of two different realiza- 
tions of quenched disorder for two sizes (a)-L = 16, (b) L=24 

of (preferably small) size L to determine the CrMES at this size. Then a new 
CrMES corresponding to a lattice of size L' > L is calculated via Eq. [2] in 
order to perform the WL simulation at this bigger size L' only in a restricted 
energy interval. 

Figure [2] shows that when the size increases new peaks which are absent for 
smaller sizes, may appear in the specific heat curve. Moreover the location of 
these peaks is strongly sample dependent. In figure [2] the specific heat curves 
of two different realisations at the same field are shown for L = 16 (a) and 
L = 24 (b). 

Figure [3] shows the energy probability distribution corresponding to the max- 
ima of Cl{T) depicted in Figures [T] and [21 It can be seen that for ho = 2.1, 
for L > 16, the probability distribution of the energy is double peaked. This 
is a signature of a first order transition, which is absent for smaller sizes. So 
the procedure proposed in [7] would be misleading if the starting sizes are 
smaller than L = 16. For the 3D RFIM L = 16 is big enough so as to render 
the use of multirange WL scheme necessary. Very recently a variation on the 
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Fig. 3. (a) Probability distribution of the energy for different lattice sizes at the 
corresponding transition temperatures T*. Double peaks appear in the probability 
distribution as L > 16 suggesting a first order transition. Probability distributions 
of different sizes have very small overlap. It shows the size of the shift from AE{L) 
to AE{L'). (b) Probability distribution of the energy for L = 24 and two different 
samples. One can get an estimate of the extended relevant energy interval for this 
size. 

WL method has been proposed that migth improve this point [17]. 

Finite size effects are well known and should be generally taken into account in 
order to obtain results for infinite systems. The additional point when consid- 
ering finite size effects in disordered systems is that changing the size implies 
a change of quenched disorder, which is obviously not the case in systems 
without disorder. This becomes particularily troublesome when large sample 
to sample fluctuations are observed. 

The observed large sample to sample fluctuations are in agreement with other 
works on disordered systems [9] [12] [18] - [20]. As the multiplicity of peaks 
of the specific heat and their locations are strongly sample dependent (see 
Figure [2]), the location of the maxima of Pt*{E), and the resulting CrMES 
depend also on the considered sample. 
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To take this aspect into account in [9] a "broadened CrMES" for WL simu- 
lation is proposed. This broadened CrMES is determined by the union of M 
CrMES each of them calculated for a different quenched random field configu- 
ration. Then, using Equation [21 the CrMES is estimated for a larger size. It is 
assumed that for the new L' > L, it will be enough to run the WL algorithm 
only in this new restricted energy interval. 

Letting aside the fact that equations [1] and [2] are not valid if the transition is 
first order, and that in general the critical exponents are not always known in 
the case of a second order transition, we will now discuss in detail the idea of 
extending to a larger size L' the CrMES calculated for L [T6]. 

As pointed out above, in the case of disordered systems the sample of size 
L' > L is a new, different sample, the corresponding Pj-. (E) may have its 
maxima located out of the CrMES that has been extrapolated via Eq. [2J In 
Figure [3]^a) one can see that there is little overlap between the probability dis- 
tributions, even for samples of slightly different sizes. In addition Figure [3](b) 
shows that this is also true for the probability distributions of two different 
samples of the same size. 

We used the DOS obtained in [lOj for the whole energy range, to recalcu- 
late thermal properties restricting the DOS to the intervals given in ref [9J. 
Let's consider first the case where a broadened interval is used. For L = 
24 the broadened interval of [9J corresponds to the energy interval 6E = 
[—2.91, —1.83]. We find that, using this broaden restricted interval, it is pos- 
sible to reproduce our results issued from the complete DOS. Unfortunately, 
though several times smaller than the whole energy range, this interval is 
located in the hardest region from the convergence point of view [10]. So 
the multirange scheme is unavoidable and several narrow energy intervals are 
needed to achieve convergence. 

In conclusion, from the efficiency and precision point of view, working with a 
very large CrMES guarantees that all the thermodynamical properties will be 
correctly reproduced, but the computational effort required is at least compa- 
rable to performing WL in the whole energy range. 

The authors of [9] are aware of this and they have proposed, as an alternative, 
to work on the CrMES of each sample. To study this proposal we use our 
data calculated in the full energy range for each sample i and we identify the 
relevant part of the energy axis for the transition of the considered sample. 
This is noted by analogy with |9j AEi. 

Conceptually AEi corresponds to the region where the energy probability 
distribution at the transition temperature is significative for the sample i and, 
in the case of a double peaked distribution, it contains the two peaks. Hence, 
we identify AEi as the energy interval where Pt*{E) is different from zero. We 
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Fig. 4. (a) Specific heat as a function of temperature T, for L = 24 , /iq = 2.1 
obtained from thermal averages using the DOS restricted to three different energy 
intervals. Interval AEij{L) corresponds to the levels given in [9] for L=24, this union 
interval reproduces a multi-peaked Cu{T) curve coincident to the one obtained with 
the whole DOS for the sample i. On the other hand the restricted intervals AEa{L) 
and AEb{L) give the curves CaiT) and Cf,(T) respectively centered around the 
temperatures Ta ~ 1.97 and T;, 2.3. These temperatures correspond to the peaks 
of the complete DOS, that gave rise to the restricted intervals, as could be expected 



then recalculate the thermal averages using only the DOS restricted to A£^j. 
For the sample shown here the Cl{T) curve we obtained using the whole DOS 
has two separated peaks. Each of them determines a T* and hence a restricted 
energy interval. 



We show that the whole structure of Cl{T) cannot be reproduced using the 
DOS restricted to the interval associated with only one of the peaks. When 
Cl{T) is re-calculated using the DOS restricted to the interval AEf (Ai?f) 
corresponding to the temperature Ta (Tb) of one of the peaks, the other is not 
reproduced (see Figure H]). 

Now let's imagine that we determine the CrMES for a lattice size L. As a 
larger size implies a new different quenched disorder, it is possible that the 
extrapolation to V > L gives a AE' which doesn't include the highest peak 
of the new sample i' of size V . In that case the results calculated for V will 
not correspond to the transition region of this new sample. 
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4 Conclusions 



We show, using tlie example of tlie 3D-RFIM, tliat the CrMES method has 
to be applied with extreme care to the study of highly disordered systems, 
because the existence of strong quenched disorder may lead to first order 
transition properties and strong sample to sample fluctuations. 

To begin with, it is worthwhile noticing that when dealing with disordered 
systems, increasing L automatically changes the distribution of the quenched 
disorder. Hence, if the interval AE{L) calculated for a sample of size L contains 
the information to reproduce the highest peak of Cl{T), due to the large 
sample to sample fluctuations, there is no warranty that the new, extrapolated 
interval for L' will reproduce the corresponding peak of C/,/ (T). Then one is 
forced to broaden the interval at L to calculate it at L'. It should be noticed, 
however, that there is no rule to control this; additionally, if the interval is 
substantially enlarged, the computational effort is similar to that required 
when one uses the WL method on the whole energy range. 

Moreover, performing a WL simulation on the whole energy range for small 
lattices, in order to determine the nature of the transition before using the 
CrMES method should be considered with some care, because it is a method- 
ology prone to reach incorrect conclusions, as we have shown in Figure [3l As 
the lattice size increases the two peaks of Pt*{E) become well separated, and 
there is a risk using this method, to seize only one of them. This migth explain 
the anomalous behaviour observed in Fig. 5 of fS] for the largest sizes studied 
in that work. 

The actual nature of the transition of the 3D-RFIM, still under discussion, 
does not change the source of the problem. Would the transition be found 
second order in the thermodynamic limit, the use of Eqs. [T] and [2] would 
be justified. Nevertheless this would neither change the problem of sample to 
sample fluctuations, nor the fact that when working with finite sizes Pt*{E) 
would still be double peaked. 



References 

[1] B. A. Berg, T. Neuhaus, Phys. Rev. Lett. 68 9 (1992) 

[2] P. M. C. de Oliveira, T. J. P. Penna, H. J. Herrmann, Braz, J. Phys 26 677 
(1997) 

[3] J.Lee, Phys. Rev. Lett. 71 211 (1993) 



10 



[4] Fugao Wang and D. P. Landau, Phys. Rev. Lett. 86 2050 (2001) 
[5] R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. 58 86 (1987) 
[6] U. Wolff, Phys. Rev. Lett. 62 361 (1989a) 

[7] A.Malakis, A. Peratzakis and N. G. Fytas, Phys.Rev.E 70 066128 (2004) 

[8] Fugao Wang and D. P. Landau, Phys. Rev. E 64 056101 (2001) 

[9] Anastasios Malakis and Nikolaos Fytas, Phys.Rev.E 73 016109 (2006) 

[10] Laura Hernandez and Horacio Ceva, Physica A387,2793 (2008) 

[11] Laura Hernandez and H. T. Diep Phys.Rev.E 55 14080 (1997) 

[12] Yong Wu and Jonathan Machta Phys. Rev. B 74 064418 (2006) 

[13] H. Behringer and M.Pleimhng, Phys. Rev. E74, 011108 (2006) 

[14] N.G. Fytas, A. Malakis and K. Eftaxias, J.Stat.Mech. P03015 (2008) 

[15] K. Binder, D.P. Landau, Phys.Rev.B 30 1477 (1984); M.S.S. Challa, D.P. 
Landau, K. Binder, Phys.Rev.B 34 1841 (1986) 

[16] A. Malakis, N.G. Fytas and P. Kalozoumis, Physica A383, 351 (2007). 

[17] A. G. Cunha-Netto, A. A. Caparical, Shan-Ho Tsai3, Ronald Dickman and D. 
P. Landau arXiv 0808-1210vl 

[18] Giorgio Parisi and Nicolas Sourlas, Phys. Rev. Lett. 89 257204, (2002) 

[19] Nicolas Sourlas, Computer Physics Communications 121-122, 183 (1999) 

[20] Shai Wiseman and Eytan Domany, Phys. Rev. Lett. 81, 22 (1998) 



11 



